例題：種子数のポアソン回帰(個体差なし)

架空植物の第i番目の個体、個体ごとに、個体サイズ\(x_{i}\)と種子数\(y_{i}\)の依存関係を見る
個体差を想定しない場合には、

glm(y~x, family = poisson, data = ...)

これによって切片と傾きの最尤推定値が得られる

GLMのベイズモデル化

サイズによる種子数のバラツキは平均\(\lambda_{i}\)のポアソン分布\(p(y_{i}\mid \lambda_{i})\)に従うとする。
この例題では個体差をそうていしていないので、ランダム効果の項はなし。

ベイズモデルの事後分布は尤度×事前分布に比例するので、以下の様に表せる。

\[p(\beta_{1},\beta_{2}\mid \textbf{Y})\propto p(\textbf{Y}\mid \beta_{1},\beta_{2})p(\beta_{1})p(\beta_{2})\]

無情報事前分布

事前に確率分布がわからない場合には、線形予測子のパラメーター\(\beta\)が\([-\infty,\infty]\)の範囲で好きな値をとって良いという無情報事前分布を設定する
無情報っぽい事前分布として扱うのは、ばらつきがかなり大きくなる様に(偏りが出ない様に)

のいずれかを使う(この本では後者を使う)

ベイズ統計モデルの事後分布の推定

winbugs使えないしstanを使う

インストール
インストールするのめっちゃめんどくせぇ。。。
インストール補足
catalinaからのインストールにはR4.0.0が最短かも
参考リンク

kubo9.stanファイルにパラメータなどベイズ統計モデルを記述
その後繰り返し回数と、burnin(最初の方の試行のうち、何回までを使わないか)を決める
なお、繰り返し回数については、複数のサンプル列(MCMCサンプル)を比較し、その大小によって収束しているかどうかを判断する

data {
  int<lower=0> N;
  real X[N];
  real MeanX;
  int<lower=0> y[N];
}
parameters {
  real beta1;
  real beta2;
}
model {
   for (i in 1:N){
      y[i] ~ poisson(exp(beta1 + beta2 * (X[i] - MeanX)));
   }
   beta1 ~ normal(0, 100); //無情報事前分布
   beta2 ~ normal(0, 100); //無情報事前分布
}
library("rstan")
## Loading required package: StanHeaders
## Loading required package: ggplot2
## rstan (Version 2.19.3, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
load("/Users/hamaokigaizaburo/Documents/kuboR/intro-statistical-modeling-master/data/d.RData")
d.dat <- list(N=dim(d)[1], X=d$x, y=d$y, MeanX=mean(d$x)) #Stanに渡すデータ
d.fit <- stan(file='kubo9.stan', data=d.dat, iter=1600, warmup=100, thin=3, chains=3)
## Warning in readLines(file, warn = TRUE): '/Users/hamaokigaizaburo/Documents/
## kuboR/kubo9_files/kubo9.stan' で不完全な最終行が見つかりました
## Trying to compile a simple C file

MCMCサンプルから事後分布を推定

d.fit
## Inference for Stan model: kubo9.
## 3 chains, each with iter=1600; warmup=100; thin=3; 
## post-warmup draws per chain=500, total post-warmup draws=1500.
## 
##         mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## beta1   1.98    0.00 0.08   1.80   1.92   1.98   2.04   2.14  1456    1
## beta2   0.08    0.00 0.07  -0.05   0.04   0.08   0.13   0.23  1535    1
## lp__  143.97    0.03 1.05 141.29 143.58 144.32 144.70 144.95  1326    1
## 
## Samples were drawn using NUTS(diag_e) at Mon May 25 16:37:16 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

これにより、各パラメータの事後分布にまつわる情報を出すことができる
ex)beta1の値の範囲は95%の事後確率で、1.8 ~ 2.12 になる
また、収束指数はRhat列、有効なサンプルサイズはn_eff列に記述されている
このとき、隣り合うMCMCステップのサンプル間の相関が高い時に、この有効なサンプルサイズは小さくなる

plot(d.fit)
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)

traceplot(d.fit)

library(ggplot2)
library(reshape2)
library(grid)

d.ext <- extract(d.fit, permuted=F)
b1.1 <- d.ext[1:500,'chain:1','beta1']
b1.2 <- d.ext[1:500,'chain:2','beta1']
b1.3 <- d.ext[1:500,'chain:3','beta1']
b1 <- data.frame(chain1=b1.1, chain2=b1.2, chain3=b1.3)
b1.melt <- melt(b1, id=c(), variable.name="chain") #ggplotで扱いやすいようにb1を再形成
b1.melt <- data.frame(b1.melt, iter=1:500)

#サンプリング過程
p <- ggplot(b1.melt, aes(x=iter, y=value, group=chain, color=chain))
p <- p + geom_line(size=0.1)
p <- p + labs(title="beta1のサンプリング過程", x="Iterations", y="")

#事後分布
g <- ggplot(b1.melt, aes(x=value))
g <- g + geom_density() + theme_bw()
g <- g + labs(title="beta1の事後分布", x="", y="")

#サンプリング過程と事後分布を並べて表示
grid.newpage()
pushViewport(viewport(layout=grid.layout(1, 2)))
print(p, vp=viewport(layout.pos.row=1, layout.pos.col=1))
print(g, vp=viewport(layout.pos.row=1, layout.pos.col=2))

複数パラメータのMCMCサンプリング

MCMCサンプリングのためのアルゴリズム

この新しい値の確率分布とは、多変量確率分布からひとつの変量をのぞいて、他の変量すべてを定数とする一変量確率分布。これを全条件つき分布という

ギブスサンプリング変化の挙動イメージ

ギブスサンプリング変化の挙動イメージ

まとめ